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Abstract. We present a systematic study of the Green functions of a one- 
dimensional gas of impenetrable anyons. We show that the one-particle density 
matrix is the determinant of a Toeplitz matrix whose large N asymptotic is given 
by the Fisher-Hartwig conjecture. We provide a careful numerical analysis of this 
determinant for general values of the anyonic parameter, showing in full details the 
crossover between bosons and fermions and the reorganization of the singularities 
of the momentum distribution function. 

We show that the one-particle density matrix satisfies a Painleve VI 
differential equation, that is then used to derive the small distance and large 
momentum expansions. We find that the first non-vanishing term in this 
expansion is always k~ 4 , that is proved to be true for all couplings in the Lieb- 
Liniger anyonic gas and that can be traced back to the presence of a delta function 
interaction in the Hamiltonian. 
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1. Introduction 

Generalized anyonic statistics, which interpolate continously between bosons and 
fermions, are considered one of the most remarkable breakthrough of modern physics. 
In fact, while in three dimensions particles can be only bosons or fermions, in lower 
dimensionality they can experience exchange properties intermediate between the two 
standard ones [1] . In two spatial dimensions, it is well known that fractional braiding 
statistics describe the elementary excitations in quantum Hall effect, motivating a large 
effort towards their complete understanding. Conversely the study of one-dimensional 
(ID) anyons is still at an embryonic stage, although recent interesting proposals for 
topological quantum computations based also on ID anyons [2] . 

In one dimension, anyonic statistics are described in terms of fields that at 
different points (xi ^ x 2 ) satisfy the commutation relations 

* A ( Xl )* A (x 2 ) = e^< x ^^ A {x 2 )^ A { Xl ) , (1) 

*a(xi)& a (x 2 ) = e- iK ^- x ^^ A (x 2 )^ A (xi) , 

where e(z) = — e(— z) = 1 for z > and e(0) = 0. k is called statistical parameter and 
equals for bosons and I for fermions. Other values of k give rise to general anyonic 
statistics "interpolating" between the two familiar ones. 

A few ID anyonic models have been introduced and investigated [3, 4, 5, 6, 7, 
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In this paper we consider the anyonic 
generalization of the Lieb-Liniger gas defined by the (second-quantized) Hamiltonian 

H= JmJ q dxd ^ A {x)d x ^ a{x) + c J dx* A (x)* A (x)y A (x)* A (x), (2) 

which describes N anyons of mass M on a ring of length L interacting through a 
local pairwise interaction of strength c (in what follow, we fix 2M = K = 1). In 
first-quantization language, the Hamiltonian is 

N f) 2 

H = -Y.q^ + 2c E S( Xi - Xj ), (3) 

i 1 l<i<j<N 

where now is the iV-anyons wave- function ^/ K (xi, x 2 , .., Xn) to exhibit a generalized 
symmetry under the exchange of particles 

* K (- • • Xj,x j+1 • • •) - e ""«(^+i^)* K (. . . x j+uXj ■ • •)• (4) 

For k = the model reduces to the bosonic Lieb-Liniger [20], while for k = 1 to free 
fermions. The physics only depends on the dimensionlcss parameter 7 = c/po (with 
Po the mean density po = N/L) and not on the two parameters separately (a part 
from obvious scaling factors). The interest in this model is mainly due to the fact 
that it is the simplest solvable model of interacting anyons and in fact it has been 
shown that it has a solution in terms of Bcthc Ansatz [3, 4]. Furthermore ID exactly 
solvable models are in a renascent age after the experimental realizations of trapped 
ID atomic gases in the last few years [21]. There is also a proposal to engineer an 
anyonic gas by trapping bosons in a rapidly rotating trap [22]. 

However, as well known for the bosonic counterpart, the Bethe Ansatz gives a 
precise characterization of the spectrum of the model and of the full thermodynamic, 
but does not allow to calculate the correlation functions. More complicated methods 
employing an algebraic formalism [23] (eventually joined to numerical calculations 
[24]) must be used in order to extract the correlation functions at arbitrary distances. 
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This with one important exception: the case of impenetrable particles, i.e. c = oo, 
that is obtainable with an anyon-fermion mapping [6] . 

Here we exploit this mapping to give a representation of the one-body density 
matrix 

P K N (x) = (^(x)nm, (5) 

in terms of the determinant of a Toeplitz matrix (completing the work started by one 
of us in Ref. [7]) whose large N asymptotic is given by the Fisher-Hartwig conjecture. 
We present a careful numerical analysis of p 1 ^ (x) and of the Fourier transform (known 
as momentum distribution function) for general values of the anyonic parameter k 
showing in full details the crossover between bosons and fermions, that is highly 
non trivial because it involves a reorganization of the singularities of the momentum 
distribution function. 

Furthermore we show that p%(x) satisfies a Painleve VI differential equation, that 
is the same as for impenetrable bosons [25] but with different boundary conditions. 
This differential equation allows a straightforward derivation of the small x expansion 
of p1f{x) that can be used to derive the large momentum expansion of the momentum 
distribution function. The first non-vanishing term in this expansion is always fc -4 a 
fact that is proved to be valid for all couplings c in the Lieb-Liniger gas and that can 
be traced back to the presence of a delta function interaction in the Hamiltonian. 

After the completion of this manuscript, a complementary approach to the same 
problem was used by Patu, Korepin and Averin [26]. They derived p%{x) in the 
thermodynamic limit as a Fredholm determinant generalizing the Lenard result for 
bosons [27]. However, despite the fact that we both consider the same correlation, the 
two approaches are complementary. In fact, our method allows to study effectively 
systems with finite numbers of particles and it is particularly suited for asymptotic 
expansions. The method of Ref. [26] allows instead for a direct generalization to finite 
temperature that in our approach is very cumbersome. The complementarity of the 
two approaches is highlighted by the fact that there is a single common formula in the 
two manuscripts (namely Eq. (43) below, derived in very different ways). We stress (as 
done in Ref. [26] ) that an explicit proof of the equivalence of the two representations 
of the anyonic correlation function is still an open problem, as in the case of bosons. 

The paper is organized as follows. In Sec. 2 we present the determinant form for 
p%(x) and derive its asymptotic behavior for large distance by means of the Fisher- 
Hartwig conjecture. In Sec. 3 we present the numerical results that allow for a 
characterization of the crossover from bosons to fermions. In Sec. 4 we prove that 
p%{x) satisfies a second order differential equation, that in the next section 5 is used 
to derive the small distance and large momentum expansions. In Sec. 6 we show that 
the power-law tail of the momentum distribution function is generally valid for the 
Lieb-Liniger model. Finally, in Sec. 7, we discuss critically our results and possible 
future investigation. 

2. Ground-state function and one-particle density matrix 

In the limit of impenetrable anyons, i.e. c — > oo, the N anyons ground-state wave 
function can be easily written down by imposing that two anyons should not occupy 
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the same position. As shown in Ref. [6], the ground-state ^ f g(x 1 , • • • , x N ) is then 



%(x u ---,x N ) 



n a & 

\<i<j<N 



Xi) 



^l(xi, ■ ■ -,x N ), 



(6) 



where ^q(xi, ■ ■ ■ , Xn) is the ground-state function of N free fermions. In the following 
we always consider the case an odd number of particles N, which corresponds to a 
non-degenerate ground state. In this case we have 

&l(xi,X2,- ■ ■ ,%n) 



and 



A( Xj 



Jtt(I-k) 



zdet 

N-l 



l.k 



N-l 



1,- 



,N, 



Xj <C X% , 

X si (Cfl • 



(7) 



(8) 



The above Anyon- Fermi mapping directly generalizes the Bose- Fermi one [28]. 
Here, the symmetry properties of the wave function are encoded in the factor 
rii<i<j<JV A(x i7 Xj) which gives the statistical phase e J7r ( 1 ~ K ) p resulting from the P 
exchanges needed for the particle positions to be brought to the ordering < X\ < 
x 2 < ■ ■ ■ < xn < L. Note that in the definition of Eq. (8), we explicitly specified 
the sign for the exchanging phase. This amounts to fix how two anyons exchange 
their positions on the ring. Moreover, in order to define the boundary conditions, we 
also fixed how a loop of one variable encircles the others. Figure (1) shows pictorially 
the example of two particles once the convention in (8) is fixed. Clearly Eq. (8) is 
equivalent to Eq. (1). 

The boundary conditions on the wave function of N anyons at positions < x\ < 
X2 < ■ ■ - xn < L can then be chosen such that (for a detailed discussion see [11]) 

• ■ ■ ,Xj + L, - ■ ■ , x N ) = e -(i-«)(^-i)-2-(i-«)0--i)^« (a;i) ..., Xj ,.. ., XN ). (9) 

A last remark is that, in general, one can allow an overall phase, coming for example 
from a non-zero magnetic flux penetrating the ring [8]. 





1 2 




2 1 



1 K K 




-m k 



Figure 1. Graphical representation of the exchange convention used in this paper 
corresponding to Eq. (8), or equivalently to Eq. (1). 



One-particle density matrix of ID anyons 5 
2.1. The one-particle density matrix as a Toeplitz determinant 

The one-particle reduced density matrix p'^ r (xi,x' 1 ) can be written in terms of the 
ground-state wave function as 

f L f L f L —K 

p%(xi,x' 1 )= dx 2 dx 3 .. dx N ^ (x 1 ,x 2 ,x 3 ..,x N )^Q(x' 1 ,x 2l x 3 , ..,x N ) . (10) 

Jo Jo Jo 

In the case of a homogeneous system Px(xi,x[) — p%(x\ — x[). We can thus set 
x[ = and study the function p%(x) with p%(0) — 1 as fixed by the normalization 
chosen in the definition (10), since the wave- function (6) is normalized to 1. From Eq. 
(9), it follows 

p%(x + L) = e^ 1 -^ N -V p %(x), (11) 

i.e. Pn( x ) is an ^-periodic function only if n = 2m/ (N — 1) with m integer, as stressed 
inRef. [7]. 

Analogously to the case of impenetrable bosons, the anyonic one-particle density 
matrix can be written as a Toeplitz determinant. Using Eq. (6), the anyonic wave 
function \I>q is 

V$(x u x 2 ,---,x N ) = -7== [I 2iA(x j -Xi)sm[TT(x j -x i )/L], (12) 

VN\L l<i<j<N 

that allows to write the integral in Eq. (10) as (in the angular variables 2nXj/L = tj) 

r 2lT r 27T 



pnV) = mj o dh '"j dtN 



N 2 

J] -A(t s - t) sin[(t a - t)/2] sin[t a /2] JJ 4sin[(tj - U)/2} 2 . (13) 

s=2 2<i<j<N 

Note that the dependence on the anyonic parameter k enters only through the function 
A(t s — t). Using the identity 

IJ 4sin[(t J -i fe )/2] 2 = [] | e ^-e^'| 2 , (14) 

2<j<k<N 2<j<k<N 

we can identify the second product in the integral (13) with the square of the absolute 
value of a Vandermonde determinant. The anyonic density matrix p 1 ^ (t) — p%(2-KxjL) 
can be finally written as 

p%(t) = Idct^ , (15) 

where the (N — 1) x (N — 1) matrix <f>£ ; has entries 

®h = l^ fj dt s e^ k - l ^^(t s ) 

= / dt s e< k - l ^4A(t s - t) sin[(t 8 - t)/2] sin[i s /2] . (16) 
^ Jo 



Note that we have 



- ^ J J dt s e^ k - l ^Asm[(t s - t)/2] sin[t a /2] , (17) 



= ^ / " A s e l(fe ^ )ts 4| sin[(t a - 0/2] sin[i s /2]| , (18) 

for /c = 1 and corresponding to the well known results for free fermions and 
impenetrable bosons respectively. 
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2. 2. The Fisher- Hartwig conjecture and the asymptotic behavior for large x 

Toeplitz matrices are fundamental objects in the study of lattice models and their 
extensive study started in the sixties for the calculation of the correlation functions 
in the classical two-dimensional Ising model. This study culminated with the Fisher- 
Hartwig conjecture [29] that relates the determinant of a Toeplitz matrix for large N to 
the analytic structure of the generating function, as briefly reviewed in the following. 
Let us consider the N x N Toeplitz matrix T with entries 

T p , q = ±J%( S y^ s ds. (19) 

(p(s) is called the generating function of the matrix. The Fisher-Hartwig conjecture is 
formulated in terms of the canonical factorization of the generating function 

R 

<p(s) = b(s) Y[ Tp r (s - s r )u ar (s - s r ) , (20) 



r=l 



where 



and 



r e -*f3(K-s+s r ) S r < S < 27T + S r , 

o< a <»„ ,21) 



u a (s) = (2- 2coss) a . (22) 

1/3 (s) takes in consideration the jump discontinuities while u a (s — s r ) encodes the 
possible singularities. The Fisher-Hartwig conjecture states that the leading term in 
the limit N ^> 1 of the determinant of T is given by 

detT ^G[b(s)] N N^r=M-^E, (23) 

where 

G[b(s)\ =e^C dsX < h ^\ (24) 

The constant E has been determined only long after the formulation of the conjecture 
[30] 

E= TT ^ e i{s r -s l )y( ar +P r )( ai -0 l ) TT G (l + a r + Pr)G{l + OV ~ gr) ^ ^ 
l<r^l<R r=l + ^ a ' r ^ 

where G(x) is the Barnes function 

G(z + l) = (27r) z / 2 e-( z +^ +1 ^ 2 )/ 2 f](l + z/fc) fc e - z + z2 /( 2fc ), (26) 

k=i 

7 B is the Eulcr constant, G(l) = G{2) = 1, and G(3/2) = 1.06922 . . .. 

The above result can be applied to the calculation of the one-particle density 
matrix of the anyonic gas, where the generating function is given by (fi K (t s ) = 
&A(t s — t) sin[(£ s — t)/2] sin[s/2] in Eq. (16). <p K {t s ) is a piecewise continuous function 
which takes the values 

= f e^ 1 "*) sin[(i s - t)/2] sin[i s /2] for < t s < t , 
[s) \ sin[(i s -t)/2]sin[i s /2] fori<i s <27r. 
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Using 



e m+i2n0 ior t s <t, 

e tf3t for t s >t, 



T (t s )T_p(t s -t)={ m r _ / (28) 



and 

2 - 2 cost, = 4sin 2 (t s /2), (29) 
the generating function can be written as 

4>-(t s ) = e- i(3t T_ K/2 (t s )T K/2 (t s - t)(2 - 2cosi s ) 1/2 (2 - 2cos(t s - t)) 1 / 2 . (30) 

Comparing Eq. (30) with Eq. (20), we have R = 2, /?i = -f3 2 = -k/2, 
a\ = a 2 = 1/2. Using the Fishcr-Hartwig conjecture we finally find 

p%{t) ~ e^ N ~ 1 ^ t N~ 1 ^- K ^ 2 [G(3/2 + k/2)G(3/2 - k/2)] 2 

X (1 - e - i *)-V4-K 2 /4-K/2(- 1 _ e i^-l/4-K 2 /4+K/2^ (3-^ 

After simple algebraic manipulations we arrive to the final result 

p K N {x) ~ (2iV)- 1 / 2 - K2 / 2 [G(3/2 + k/2)G(3/2 - k/2)] 2 

-1/2-k 2 /2 



x e i7r K (AT :C /L-l/2) 



sm( T J , (32) 



where we have reintroduced the variable x via t — 2ttx/L. 

For k = 0, Eq. (32) reproduces the famous Lenard result [31] for impenetrable 
bosons. For k = 1 instead it does not reduce to the free fermion result p\[{x) = 
sm(NTTx/L)/(Nsm(TTx/L)). In fact, as shown, in Ref. [10], Eq. (32) can only get 
the mode at k = —hp = —{N — l)/2, but not the other one at k = kp, that only for 
k = 1 is degenerate with the other. This is not inconsistent, because for k = 1 the 
Fishcr-Hartwig conjecture does not hold. 

All the results up to this point correspond mainly to an extended and detailed 
version of those already appeared in the letter [7] . All new results are reported in the 
following. 



3. Explicit calculation of the reduced density matrix and momentum 
distribution function 

Eq. (15) provides a representation of p%{x) that can be easily used to derive the 
correlation function at finite N. In particular, for a small number of particles 
we have explicit expressions of p^(x). For N = 3, using the variable t = 2irx/L, 
< t < 2?r, wc find 

PN= 3 (t) = ^ [Stt 2 - & (-15& + 87rf - 2£ K i 2 ) + & cos2t 

+4(27r-^(t + 2))(27r-^(i-2))cost+12^(27r-^)sint] , (33) 
where £ K is defined by 

& = 1 - e l{1 - K)7r . (34) 

In the same manner and with the help of Mathematica, we can expand the determinant 
up to N — 11. However, the formulas for general k are too long to be reported here. 
Larger values of N can be easily worked out numerically The small N analytic 
expressions are then practical cases to check numerical calculations, asymptotic 




Figure 2. Momentum distribution function ^^(fc) (i.e. N = L = 81) for 
several statistical parameters re showing how the peak decreases with increasing 
re and becomes the Fermi discontinuity at k = —kp. Inset: Zoom in the region 
kp < k < 2kp showing that the weak singularity at k = 2kp for re = becomes 
the discontinuity of the function at k = kp for re = 1. 



expansions etc. as we will extensively do in the following. For example, by fixing 
t = 2tt in Eq. (33), and in the analogous ones not reported here, we can explicitely 
verify that p% =3 (L) = (1 - Z K fp% =3 (0) and p% =5 (L) = (1 - £ re ) 4 p& =5 (0) etc. in 
agreement with Eq. (11). 

We now focus on the properties that can be derived from the numerical 
calculations of p%{x). The correlation function in the real space has been already 
worked out numerically for several k in Ref. [7]. The agreement with the Fisher- 
Hartwig result Eq. (32) has been shown to be excellent for small k. When k get 
closer to 1, by using an harmonic fluid approach, it has been shown [10] that a second 
mode not captured by the Fisher-Hartwig, becomes important and is responsible for 
the right behavior at the fermionic point (these results were later rederived with a 
different approach [11]). Overall we can safely state that the main qualitative features 
of p%(x) that can be extracted by numerics are well understood (with the important 
exception of fixing the amplitude of the second mode, see the discussion section for 
details). The same is not true for the momentum distribution function n^(k) which 
has only been considered marginally [7, 10]. 
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-1 1 

k/(N-l) 

Figure 3. Scaling of the momentum distribution function n l ^°- 1 (k) with varying 
N = 21,41,61,81,101,121. Upper inset: n K N {k) as function of k. Main Plot: 
Scaling of the peak by plotting n^ r (fc)_/V Q " with a K = 1/2+k 2 /2 versus k/(N— 1). 
Central inset: n^(k)N versus k/(N — 1) that gives the collapse for large k, and 
works well also for moderate values of k. 



Here we fill this gap. For practical purposes we only consider the case of a periodic 
p%{x), i.e. K — 2m/ (N — 1), so that n%(k) can be defined as the Fourier transform 

n K N {k) = ±J dxe 2 " ik */ L p%(x), (35) 

with k integer. We calculated nfj(k) for all the values of n £ [0, 1] that give an L- 
periodic Ppj( x ) for N = L = 21, 41, 61, 81, 101, 121 (higher values of N could have been 
easily obtained, but the values considered are sufficient for our aims). The method is 
very simple: we calculated numerically p%(x) for enough equispaced x's and then we 
considered the fast Fourier transform of these data, obtaining very accurate results. 

Let us first discuss the more natural question we can address with this method: 
how with changing n we can smoothly interpolate between the impenetrable boson and 
the free fermion momentum distribution functions. We remember to the reader that 
these two correlation functions are very different. The free fermion one is obviously a 
Fermi-Dirac distribution 

w M fl/N for|fc|<fc F , 
HN{k)= 1 for|fc|>fc F , (36) 
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k/(N-l) k/(N-l) 

Figure 4. Scaling of the momentum distribution function n^ 5 (fc) with varying 
N = 21,41,61,81,101,121. Upper inset: n K N {k) as function of k. Left Plot: 
Scaling of the peak by plotting n^ r (fc)iV Q ' t with a K = 1/2 + /t 2 /2 versus k/(N-l). 
Right Plot: jijy(fc)iV versus k/(N — 1) that gives the collapse for large k. 

with kp = (N — l)/2. Oppositely at k = 0, n^(fe) has the characteristic of a strongly 
interacting system with a peak in 0, going like TV -1 / 2 [27] and a large-momentum 
power-law tail [32, 33]. Both the limits are symmetric for k — > — fc. 

Fig. 2 shows Wgi(fc) (we take iV = 81 as a typical large enough value) for several 
k between and 1 (the figure with all the possible k making p%(x) periodic has 
too many curves to be readable). The peak at k = for n — shifts backward at 
k = -k(N- l)/2 and its height decreases from N^ 1 / 2 to N~ a * with a k = 1/2 + k 2 /2. 
Both these features are encoded in the Fisher-Hartwig result Eq. (32). Approaching 
k = 1 this peak becomes the discontinuity at k = —kp. However, the Fisher-Hartwig 
result only applies very close to this maximum value and cannot explain how the 
discontinuity at k — kp is produced, that instead is well understood from numerics. 
For bosons, it is known from the pioneering paper of Vaidya and Tracy [34] that not 
only k = is singular, but there are additional weaker singularities at all the points 
k = 2mkp with m integer. For example, at k = ±2kp the second derivative of n ( j v (k) 
is divergent in the thermodynamic limit, but this singularity is so weak that is hardly 
seen in Fig. 2. Increasing k, all these singularities move backward of —k(N— l)/2. In 
particular, the first one at 2kp moves at 2kp — n(N — l)/2 and becomes sharper. This 
is emphasized by the inset of Fig. 2, where we zoom close to this region showing how 
the derivative of n^(k) develops a larger discontinuity increasing re, that becomes a 
discontinuity of the function itself for k = 1 at k = kp . 
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-1 -0.5 0.5 1 

k/(N-l) 

Figure 5. Scaling of the momentum distribution function n^ 9 (fc) with varying 
N = 21,41,61,81,101,121. Upper inset: n^(fc) as function of k. Main Plot: 
Scaling of the peak by plotting n^ r (fc)A rQ " with a K = 1/2 + k 2 /2 versus k/(N — 1). 
Central inset: n l ^ f (k)N versus k/(N — 1) that gives the collapse for large k. 

We find that this mechanism to connect smoothly n = to n = 1 is very 
interesting because of its simplicity. Furthermore we believe that it should be valid for 
the Lieb-Liniger model at arbitrary coupling, but unfortunately it is still impossible 
to have the analytic structure of the singularities in the general case. However, this 
structure is compatible with the results from the harmonic fluid approach [10, 11], 
where the amplitudes that fix the strength of the singularities are free parameters. 
Finally it is likely that this mechanism would be valid also for other models of 
interacting anyons. 

There are other interesting features that can be extracted from our numerical 
calculations. We notice that, in contrast with the cases k = 0, 1 the momentum 
distribution function is highly asymmetric. Increasing k the decay to the left of the 
peak becomes more rapid, while for k > — n(N— 1)/2, n^(fc) slowly develops a plateaux 
that becomes the Fermi sea at k = 1. However for value of \k\ larger than kp the 
momentum distribution function tends to restore the symmetry k — > —k. This is clear 
from Figs. 3, 4, and 5 where we plot n^(k) for k = 0.1,0.5,0.9 in logarithmic scale 
to magnify the small values the functions take. In the following sections we will show 
rigorously that these tails follow a k~ A power-law and are symmetric with respect to 
k = (the first asymmetric term is only at level k~ 9 when k — 2m/ (N — 1)). 
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The three Figs. 3, 4, and 5 show also other interesting features of the momentum 
distribution function. We plot n^(k)N aK with a K = 1/2 + k 2 /2 as function of 
k/(N — 1) = 2k I k F that should show a perfect data collapse for all the values of 
k well described by the Fourier transform of the asymptotic Fisher-Hartwig result 
valid for large x, and so close to the peak. It is evident from the figures, that the 
collapse is effective only extremely close to the peak. Oppositely for large k, the 
momentum distribution function is not expected to show any anomalous scaling and 
to be proportional to l/N (this will be shown rigorously for large k in the following 
sectionsL note the different normalizations of n^(x)). For this reason we also plot 
Nn^(k) that in fact shows a good collapse of the data for large k. But not only: also 
for intermediate values of k the various curves fall on the same "scaling" function. 

4. Differential equation for p l j i (x) 

An effective way of characterizing correlation functions for ID strongly interacting 
systems is to find a differential equation that the correlation satisfies and from this 
extract the analytical properties of the solution. For the ID impenetrable Bose gas 
this program started with the work of Jimbo et al. [35] proving that p\ {x) (i.e. in 
the thermodynamic limit) satisfies a second order Painleve differential equation of the 
V kind Py . This result allowed to derive several terms in the asymptotic expansion 
for large distances [35]. Later Forrester and collaborators were able to show that also 
the finite N correlation function p%(x) satisfies a Pyi Painleve equation and they 
pointed out useful connections with the theory of random matrices. Exploiting this 
connection, we show here that also the anyonic one-particle density matrix pfj (x) can 
be characterized via a second order non-linear differential equation Pyi- 

The starting point is that \^q(xi,X2, ■ ■ ■ xn)\ 2 can be viewed as a probability 
distribution function for some class of random matrices [36] (^J(xi, x 2 , ■ ■ ■ xn) is the 
free fermion ground state). For periodic boundary condition, the appropriate matrix 
ensemble is the so called unitary circular ensemble [37], for which 

[^(x!,^,...,^)! 2 ^^^ [] 4sin[^(x J - a;fc )/L] 2 -Ev(U(N)), (37) 

l<j<k<N 

where Ev(U(N)) is the eigenvalue probability distribution function of unitary matrix 
U(N) with uniform measure. Eq. (13) for p%(x) can thus be interpreted as an average 
in the circular unitary ensemble 

N 2 

P%(t) = (II ~ A ( ts - *) sin ^ - *)/ 2 l sin[t s /2]) Ev(u(N)) . (38) 

It is known (see for instance Ref. [38] and references therein) that the averages in 
various random matrix ensembles are related to the Painleve differential equations 
which are classified in six types: Pj,---Pyi- In the case of the average (38) it 
has been shown by Forrester et al. [25] that one can evaluate p^ via non-linear 
differential equation of Pyi type. In fact, P2{x), P§(x), • • • Pn+i( x ) can be considered 
as a sequence of functions T3[1](u),t 3 [2](u), ■ ■ ■ t 3 [N](u) (we introduced lnu = 2ix) 
which represents one of the so called r-function sequence occurring in the Pyi 
systems. The first function of this sequence, T3[l](u), turns out to satisfy a Gauss 
hypergeometric equation, which admits two independent solutions. The parameters 
appearing in the hypergeometric equation do not depend on the anyonic parameter k. 
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The effect of the statistics is to select a particular solution in the bidimensional space 
of solutions of the hypergeometric equation. In other words, the anyonic statistics 
enters only in the determination of the boundary conditions. Using the so called 
Backlund transformations, which leave the form of the Pyi equations unchanged, 
one can systematically construct all the T3 [N] function from the T3 [1] and verify that 
the form of the corresponding differential equations do not depend on the statistics. 
Specifically, the function (J%{u) (directly related to t^[N]{u)) 

lnu 



o-%{u) = u(u - l)]np%(x(u)), 



x(u) = 



2i ' 



satisfies the following second order differential equation 



u 2 (u-l) 2 



r 



du 2 



o- N (u) -(u - l)-^a%(u) + 1 



(39) 



(40) 



du 



a K N {u) 



N 2 



a K N {u)-{u-l)—a K N {u) 



= 0, 



that in fact does not depend explicitely on the anyonic parameter k. The study of 
the boundary conditions for general values of k is missing in Ref. [25] , where only the 
bosonic and fermionic cases (k = 0, 1) have been discussed. The general k dependence 
of the boundary conditions (£ K is defined in Eq. (34)) 



lim <jn{u) 



U->1 



N 2 -l 
12 



(u-lf 



N 2 



1 



24 



N{N 2 - 1) 
48^ 



(u If 



(41) 



are derived by expanding p^{t = 2nx/L) in power of x. Let us show how to derive 
this result. Using Eqs. (8) and (34), Eq. (13) can be written as 

1 N ( f 27r /"* \ N 2 

p%(t) = ^ n {j dts - & j dts ) n ^ sin[(<s - i)/2] sin ^/ 2 i 

x [] 4sin[(t,-t l )/2] 2 . (42) 

2<i<j<N 

Expanding the above expression in powers of £ K , p%(t) is expressed in terms of the 
(n + l)-particle density matrices p} N ^ (t, t 2 , ■ ■ ■ t n+ \\ 0, ti, • • • t n +i) of a system of N 
free fermions 



^ — —. — / dt 2 dt 3 --- dt n+ \p} nN s{t, t2, ■ ■ ■ tn+i; 0, t2, ■ ■ ■ i n +i 
n=1 n - Jo Jo Jo 



), (43) 



The n— th term of the above expansion is proportional to the n— th power of t. Using 
the Wick theorem, the (n+ l)-particles density matrix p]^(t, t 2 , ■ ■ ■ t n+ i; 0, t2, ■ ■ ■ t n +i) 
is expressed as a products of one-particle density matrices p\,. Using 

sm(Nirx / L) 



pl ^ x) ~ Nsmfrx/L)' 
the 2 free fermions density matrix pl 2 N \(t, t 2 ; 0, t 2 ) is 

p\ 2 N) {x,X2]0,x 2 ) = p 1 N (x)p 1 N (0) - p 1 N (x 2 )p 1 N (x - x 2 ) = 

s'm(Nirx/L) s'm(Nirx 2 /L) sin( N it (x — x 2 )/ L) 
N sin7rx/ L N sm(irX2 / L) N sm(ir(x — x 2 )/L) 



(44) 



(45) 
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At the first order in £ K , we have 
sin(7V7rx/L) 



Pn{x) = 



N sui(ttx / L) 



_k f 

L Jo 
which gives 



dx 2 N 



sm{-Kx/L) 



sm{Nirx/ L) s\n{Nitx 2 / L) sin(A^7r(a; — x 2 )/L) 
sin(7rx 2 /L) sin(7r(x — x 2 )/L) 

1) 



0(^),(46) 



lim p%(x) = 1 



TV 2 



-(ttx/L) 2 + £ k 



18tt 



(ttx/LY 



(47) 



that is equivalent to Eq. (41). 

Note that Eq. (42) shows explicitely that the non real (i.e. the imaginary parts) 
terms in p%(x) come from £ K alone. This is suggestive for an artificial interpolation 
between bosons and fermions that does not require complex phases. 



5. Short distance and large momentum expansions 

Starting from the work of Jimbo et al. [35] differential equations satisfied by correlation 
functions become one of the most powerful tool to obtain asymptotic expansions. Eq. 
(40) joined with the boundary condition (41) allows to obtain the small x expansion of 
p%(x), on the same line of Ref. [39] for impenetrable bosons. In fact, by substituting 
a small x power series for p%{x) into the differential equation, we obtain equations 
which define all but one of the coefficients. In particular the resulting equation for the 
coefficient x 3 vanishes identically, and to fix this parameter we require the boundary 
condition (41). With the help of Mathcmatica, we found straightforward to obtain the 
first 25 terms of p%(x), but it would require far too much space to exhibit all these 
here. In the variable t = 2nx/L up to order i 10 we find 

^ = 1 -V + ^V' (48) 

(3TV 2 -7)(7V 2 - 1) 4 iV(lliV 2 -29)(7V 2 -l) 5 
+ 2*360 1 2 5 2700tt t 

(37V 4 - 187V 2 + 31)(iV 2 - 1) 6 7V(183iV 4 - 12107V 2 + 2227)(7V 2 - 1) 7 
2615120 1 2 7 1587600tt 1 

(TV 2 - 1)[(157V 6 - 1657V 4 + 7177V 2 - 1143)tt 2 + 567V 2 (7V 2 - 4)£ 2 ] 8 
+ 2 8 5443200tt 2 1 

7V(-22863+ 138677V 2 - 30177V 4 + 2537V 6 ) (TV 2 - 1) 9 
2 9 142884000tt 1 

TV 2 - 1 

tt. ^[525(37V 8 - 527V 6 + 4107V 4 - 16367V 2 + 2555W 2 + 

2 10 62868960000tt 2 L 1 ' 

887V 2 (4897V 6 - 46067V 4 + 127617V 2 - 8644)£ 2 )t 10 . 

In the bosonic case, k — (£ K = 2), we find the result obtained in Ref. [39]. In the 
fermionic = 1 (£ K = 0), the above formula gives the small x expansion of p 1 N (x) 

ofEq. (44). 

Note that the dependence of the previous expansion up to the order t 7 is trivial, 
in the sense that the even terms are the same as for bosons, while the odd ones only 
get a global factor £ K . While the latter property is true at all the orders we calculated, 
the former has a non-trivial £ K dependence that firstly shows up in the term t 8 with 
a factor proportional to £ 2 . Higher powers in t show also higher powers of £ 2 . 
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5.1. Large moment expansion 

Because of the periodicity properties of p%(x) in Eq. (11), the momentum distribution 
function cannot be defined simply as a Fourier series. Its definition should be changed 
according to (for convenience we also introduced a factor 1/N compared to the 
definition in Eq. (35)) 



Pn( x ) = x ^2 n K N (k)exp\-2ni(k + S K )j], (49) 

n— — 

where the shift S K is defined by 



-oo 



5* = {l + (K-l)^}, (50) 

and {x} = x— [x] stands for the non-integer part of x. This shift in the definition does 
not matter when p%{x) is periodic, i.e. for n — 2m/(N — 1), when 5 K = 0. Inverting 
Eq. (49), the momentum distribution function is 

N f L r xi 

n%(k) = — J cfoexp \2m(k + 6 K )-\ p%(x), (51) 

that gives the probability occupation of the states with momentum 2ir(k + 5 K )/ L. For 
a small number of particles it is possible to obtain close expressions of n^(fc). For 
N = 3, from Eq. (33) we have (we stress that this formula is valid only for < n < 1) 

nK , k , = 4 cos 2 (ttk/2) (4 sin (ttk) p 2 (k, k) + Trp 7 (k, k)) 

7T 3 (-4 + fc 2 +2fcK + K 2 )(-fc + fc 3 -K + 3fc 2 K + 3fcK 2 + K 3 ) 3 ' 

where the p2{k,K) and p 7 {k,K) are respectively polynomials of order 2 and 7 in k 
p 2 (k, k) = 7k 2 + Ukn - 1 + 7k 2 , 

p 7 (k, k) = 3fc 7 + 21fc 6 K + 7fc 5 (9K 2 - 2) + 35fc 4 K(3K 2 - 2) + 7fc 3 (l - 20k 2 + 15k 4 ) + 
+7fc 2 K(3 - 20k 2 + 9k 4 ) + fc(4 + 21k 2 - 70k 4 + 21k 6 ) + k(4 + 7k 2 - 14k 4 + 3k 6 ). 
For k = 0, this reduces to the bosonic distribution function n§(fc) [39] 



"3(0)4 + ^ ng(±l) = |, ng(±2) 



35 



3 2tt 2 ' , "3V^ i /- 3 ' 3V ;_ 367r 2 ' 
_ 2(3F-14fc^ + 7fc 3 + 4fc) 
"3(*J- (_4 + P)(-A : + fc 3 )% 2 ' |fc|>2 ' (53) 
and for k = 1 to the Fermi distribution 

n 1 M=n\(±l) = l, nl(k)=0 \k\ > 1 . (54) 
The close expression for nj?(fc) and n^{k) are too long to be reported here. 

5.2. Large momentum asymptotic expansions 

From the exact solutions for different values of N (N = 3, 5, 7) and k, we can study 
the behavior of the distribution of large momenta k ^> 1 by studying the asymptotic 
of the function n^(k — S K ). For N — 3 we find 

6cos 2 (ttk/2) 14cos 2 (ttk/2) 22 cos 2 (ttk/2) 56 cos 2 (ttk/2) sin(7TK) 

30cos 2 (W2) , 384cos 2 (W2)sin(^K) 12 
+ ^27To + ^37Ti + )• I 55 ) 
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Note that, for generic values of k, we have odd terms appearing in the asymptotic 
expression (55). As another illustrative example of the large k behavior of n^(fc), we 
give the large k expansions for n®(k), n^ 2 (k) and n^ 3 (k) 
50 410 2570 14234 ,,_ 12 . 

n 5y k > - ^4 + ^2js + ^8 + ~ k ~W~ + °( k )' 

i/2 25 205 1285 4900 7117 150960 n ,,_i 2 x 

i/3,, , /on 75 615 3855 3675\/3 21351 n . 

For N = 7 the large k expansions for different values of k (giving periodic Pn( x )) 
take the form 
196 3332 44604 540820 n( ,,_i 2 , 

n7[k) ~^¥ + ^ + ^W + ^k^ + 0[k h 

!/3/ln _ J£L 33453 86436^3 405615 5768280^3 _ 12 

2/3 _ 49 833 11151 28812^3 135205 1922760^3 , _ 12 . 
717 () ~^ + ^ + ^fc^ + 7r3fc9 + ^k^ + ^3jfcU )• 

The large expansion of n^(fc) for general values of k and TV can be obtained by 
means of the Mcllin transform. Given a function /(x), the Mellin transform is defined 

as 

<P(s) = / dxx^fix). (57) 
Jo 

whose inverse is 



1 



c+oo 



/(*) = ^- / dsx- s 0( S ), (58) 



c— zoo 



where the above notation implies a contour integral taken over a vertical axis in the 
complex plane in the corresponding fundamental strip. 

The Mellin transform establishes a direct mapping between the asymptotic 
expansion of a function f(x) near x = and the set of singularities of the transform 
4>(s) in the complex plane. This technique has been already applied in Ref. [39] to 
the bosonic case (i.e. k = 0). 

The asymptotic expansions we compute from the small N exact solutions suggest 
the following general large k expansion for n^(fc — S K ) 

K /, ex a l b l a 2 b 2 ft 3 ^3 

n ^ k -^^^ + ^ + ^ + i + i + ^ + -" (59) 

We use the Mellin transforms of the cos(27rnx) and sin(27rnx) functions 
cos(2nkx) = — I dsT(s) cos(tts /2)(2nnxy s < c < 1, 

Z7TJ J c—ioc 

sin(27rfcx) = — / dsT(s)sm(irs/2){2irnxy s - 1< c < 1. (60) 
Ziri Jc-iso 

Plugging Eq. (59) and Eq. (60) in Eq. (49), we have 

1 1 r c + io ° 

Re [p%] (x) = -n K N (0) + — ds F(s) cos(irs/2)(2Trx)- s g(s), 

L 2m J c _ ioo 

~\ rc+icc 

Im [p%] (x) = — dsr(s) S m(Trs/2)(2irx)- s g( S ), (61) 

2m Jc-ioc 
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where 

oo oo 

g(s) = geven(s) + g odd ( S ) = £ a^(2j + 2 + s) + ]T &K( 2 ^ + 3 + s )> ( 62 ) 

j=i i=i 

with r(x) and Ci x ) the Gamma and Ricmann function respectively. By closing the 
contour of (61) on the left, the above integral is given by the sum of the residues of 
the functions T(s) cos(7rs/2)g(s) and T(s) sm(n s / 2) g(s). The poles of these functions 
are all simple and located at the points s = 0, —2, —3, —4, —5, • • • for the integral with 
the cosine and s = —2, —3, — 4 • • • for the integral with the sine. 

In the integral (61) the singularities arise for s — —2m, m > from the 
r(s) cos(7rs/2) function and from the m — 1-th term of the sum g odd (s) while f° r 
s = —(2m + 1), m > 1, from the Riemann function appearing in m-th term of the 
sum (feven(s). Defining the function f c (z) as 

/ c (-2m) = Rcs[r(s) cos(7rs/2).g(s), s = -2m], m > 1, (63) 

the resulting expansion then reads (t = 2-kx/L) 

Be\o K }(x) - 1 - ht2± t 2 + ^t 3 + /c( ~ 4 V - ^ f 5 ^ 

Re[p N \(x)-l m t + m t + im t i2QN t (64) 

/c(-6) 6 ™g ,7 /c(-8) 8 ™l 9 , O / f 10N f65) 

64iV 5040iV 2567V 3628807V 1 '' K ' 
The coefficients of the even terms in the above expansion contain then all the terms 
of the sums in g(s). Note that the term arising from the residue at s = combines 
with n^(0)/X to give p%(0) = 1. The first terms for example are 

/ c (-2) - 4 5 (-2), 

/c(-4) = ^ ^even(-4) + £ &«<(2j + 3 + s) j + (66) 

Comparing Eq. (48) with Eq. (65), we have the coefficients 
N 2 (N 2 - 1) 

N 2 (N 2 -l)(-29+UN 2 ) 
° 2 " 720^ RC[U 

K N 2 (N 2 - 1)(2227- 12107V 2 + 183JV 4 ) 
° 3 = 40320^ Rc ^]> 

K TV 2 (JV 2 - 1) (-22863+ 138677V 2 - 30177V 4 + 2537V 6 )) 

° 4 201600^ Re[ ^ ] - (67) 

For TV — 3, 5, 7 we recover the exact results given above. 

In an analogous way the coefficients are derived from the sine integral. Here 
the singularities for s = —2m, m > 2, arise from the Ricmann function appearing in 
(m — l)-th term of the sum g odd {s) while the singularities for s = —(2m + 1), m > 1, 
arise from the T(s) s'm(ws/2) function and from the m-th term of the sum <?even(s). 
Defining the function f s (z) as 

/ s (-2m-l) = Res[r(s)sin(7rs/2) 5 (s),s = -2m-l], m > 1,(68) 
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we find the following expansion (t = 2nx/L) 

Im^j^j - ^ t + ^t + ^ N t 72QN t + 

+ 1287V + 201607V + 5127V 9072007V + 1 j ' 1 j 
Finally, comparing Eq. (48) with Eq. (69), we have the coefficients &^ 

bl = bl = 0, 

567V 3 (7V 2 -1)(4-57V 2 +7V 4 ) 2 
• 3 ~ 34560^ lm ^' 

K = _ 7V 3 (7V 2 - 1) 2 (8644 - 41177V 2 + 4897V 4 )) 2 
4 201600tt 3 ^ 1 ; 

We have shown that the first odd term appearing in the large k expansion of n^(fc) 
is the one at the order fc~ 9 . Again, one can verify the above values match with the 
exact results for a small number of particles (TV = 5, 7). These odd terms are very 
important because they represent the onset of the asymmetry of n^(k) for large k. 
Note that the first non-vanishing odd term is k~ 9 because in the expansion (48) the 
first non trivial terms shows up at i 8 . 



6. Large k asymptotics for a finite interaction anyonic gas 

In the case of the Bose gas, it is well understood [33] that the k~ 4 tail of the momentum 
distribution function (and the corresponding small x behavior) are not a prerogative 
of the impenetrable limit, but are a signature of the delta- function interaction and 
so are general features of the Lieb-Liniger model, as confirmed by direct numerical 
calculations [24]. It is easy to generalize this result to the anyonic case and find the 
same result. 

In order to prove this, we need to briefly introduce the Bethe Ansatz solution 
of the Lieb-Liniger gas [3, 4]. For finite arbitrary coupling c, the eigenstates 
%?%(xi, . . . , xjy) can be written as [3] 

^(x u ...,x N ) = <S> K (x ll ... 1 x N )^° £ (x ll ... 1 x N ), (71) 

where the phase function 

$ K (Z1, ...,X N ) = e -^N(N-l)n/4 -Q e «r/«(*,-x p )/2 _ ^ 

q>p 

encodes all the statistic dependence of the wave functions (we adapted the anyon 
convention used in Ref. [3] to ours). In this way the problem is equivalent to find the 
bosonic wave-function ^-(x^ . . . ,xn) (in fact this method is usually called Anyon- 
Boson mapping, oppositely to the Anyon-Fermion mapping exploited in the rest of 
the work and valid only for c = oo). 

Kundu [3] showed that the bosonic coupling c is c = c/cos(7tk/2) and so all 
the thermodynamic quantities only depends on this effective bosonic coupling. For 
example, limiting to the L-periodic cases, we can write the ground-state energy per 
particle as 

^=p 2 e( 7 , K ) = p 2 (7, K = 0), (73) 
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where p a — N/L is the mean-density, 7 = c/ p , analogously 7, and e(7, k) is implicitly 
defined above. For k = 0, e(7, 0) has not a close analytical expression, but it is well- 
known and tabulated [20]. 

In Ref. [33] the large k behavior of the momentum distribution function 
for the bosonic Licb-Linigcr model was simply derived by the fact that the wave 
function at the point of contact of any two particles undergo a kink in the derivative 
proportional to the value of the cigcnfunction itself. The modification in the anyonic 
case is a straightforward consequence of Eq. (71): only the imaginary part of the 
eigenfunction %!"(xi, . . . , Xn) is discontinuous at a point of contact (and being odd 
does not contribute to the leading term for large k, in analogy with the well-known 
fcrmionic case), while the real part is continuous and its derivative has the kink of 
the corresponding boson wave-function ^(xi, . . . ,xn). Consequently the coefficient 
of the large momentum tail is the same as the one of the bosonic model at coupling 



where we adapt the result of Ref. [33] to the normalization p 1 ^ (0) = 1 . Considering 
the large 7 expansion e'(j, k) — cos(ttk/2)Att 2 /(3j 2 ) [4] we agree with the result of 
the previous section in the impenetrable limit. 

7. Discussions 

In this paper we presented a systematic study of the momentum distribution function 
and of the one-particle reduced density matrix of the anyonic generalization of the 
Lieb-Liniger model, obtaining the first full analytic description of the crossover from 
bosons to fermions for a strongly interacting model of anyons. Particular attention 
has been devoted to the large momentum and small distance behavior that can 
be analytical obtained with a proper generalization of the methods employed for 
impenetrable bosons. In the complementary regime (large distance) we obtained the 
leading term by applying the Fisher-Hartwig conjecture. Corrections to the leading 
behavior have been investigated numerically. In this regime, we find the evidence 
that the second mode of the harmonic fluid approach [10] plays a fundamental role 
in describing the correct crossover to the free fermionic regime close to the point 
k = kp. This calls for further analytical studies of the singularities of the momentum 
distribution function. Several techniques could be used to tackle this problem. In 
the bosonic case, the systematic large x expansion has been mainly exploited via the 
solution of the second order differential equation Py [35, 39], but this required the 
knowledge of the proper boundary condition for large x that was available only from 
the mapping to the lattice XX model [34] whose equivalent is not yet known for anyons. 
An alternative method, that has been applied successfully to bosons [40], would be 
to consider the replica approach. Work in this direction are in progress. Another 
effective way to obtain the asymptotics (maybe even beyond the impenetrable limit) 
would be the approach of Ref. [41]. 

The knowledge of the anyonic correlation functions beyond the impenetrable limit 
is instead still very limited. Only the power-law structure (and the corresponding 



c, i.e. [33] 




(74) 
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singularities) are known from the harmonic fluid approach [10] (or equivalently from 
conformal field theory [11]). A general approach would be to generalize the quantum 
inverse scattering methods for bosons [23] to anyons and then to mix integrability and 
numerics to get the full correlation functions from the form-factors (on the line of Rcf. 
[24]). But this is a very ambitious and long project. However, there is an interesting 
regime where this is maybe not necessary: in Ref. [4], exploiting the correspondence 
c = c/ cos(7tk/2) [3], it has been argued that for 1 < k < 2 the corresponding bosonic 
model is attractive, independently of the coupling constant c. The attractive Bosc 
gas is known to form bound-states whose wave-functions are known with exponential 
precision [42] . This allows the explicit analytic calculation of the correlation functions 
[43], that can eventually have some interpretation in anyonic language. 
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